rm(list = ls())


# Parameter
W1.F <- 1
W1.I <- 0.7
V.bar <- 2.25
Lambda <- 0.6
Cost <- 0.5


# Functions
w.0 <- function(v){
  out <- exp(-v)
  return(out)
}

LHS <- function(v, v.bar = V.bar){
  q <- (1 - pnorm(v.bar, mean = 1)) / (1 - pnorm(v, mean = 1))
  out <- pbeta(q, shape1 = 2, shape2 = 2)
  return(out)
}

RHS <- function(v, 
                cost = Cost, lambda = Lambda,
                w1.F = W1.F, w1.I = W1.I){
  top <- cost - lambda*(w1.I-w.0(v = v))
  bottom <- w1.F - w1.I + lambda*(w1.I-w.0(v = v))
  out <- top / bottom
  return(out)
}

find_equilibrium <- function(cost = Cost, v.bar = V.bar,
                             lambda = Lambda,
                             w1.F = W1.F, w1.I = W1.I){
  v.demand <- -log(w1.F - cost)
  f <- function(x){LHS(v = x, v.bar) - RHS(v = x, cost = cost, lambda = lambda,
                                           w1.F = w1.F, w1.I = w1.I)}
  return(uniroot(f, lower = v.demand, upper = V.bar)$root)
}

# Figure
V <- seq(0, 5, length.out = 1000)

pdf('./figures/fig_d2.pdf', width = 14, height = 7)
par(mfrow = c(1, 2))
plot(V, LHS(v = V), 
     type = 'l', lwd = 5, col = 'gray', 
     main = 'Informal Wage High', cex.main = 1.5, 
     xlab = '', ylab = '', xaxt = 'n')
title(xlab = 'Application Threshold', cex.lab = 1.5, 
      line = 1.25)
title(ylab = 'Probability / Adj Costs', cex.lab = 1.5, 
      line = 2.25)
lines(V, RHS(v = V), lwd = 5)
abline(v = -log(W1.I), col = '#666666', 
       lty = 3, lwd = 4)
segments(x0 = find_equilibrium(), y0 = 0,
         x1 = find_equilibrium(), 
         y1 = LHS(v = find_equilibrium()),
         lwd = 4, lty = 2)
legend(x = 2.3, y = 0.8, lwd = c(3, 3), col = c('gray', 'black'),
       legend = c('Admission Probability', 'Costs / Stakes'), 
       bty = 'n', cex = 1.3)

plot(V, LHS(v = V), 
     type = 'l', lwd = 5, col = 'gray',
     main = 'Informal Wage Low', cex.main = 1.5, 
     xlab = '', ylab = '', xaxt = 'n')
title(xlab = 'Application Threshold', cex.lab = 1.5, 
      line = 1.25)
title(ylab = 'Probability / Adj Costs', cex.lab = 1.5, 
      line = 2.25)
lines(V, RHS(v = V, w1.I = 0.1), 
      lwd = 5)
abline(v = -log(0.1), col = '#666666', 
       lty = 3, lwd = 4)
segments(x0 = find_equilibrium(w1.I = 0.1), y0 = 0,
         x1 = find_equilibrium(w1.I = 0.1), 
         y1 = LHS(v = find_equilibrium(w1.I = 0.1)),
         lwd = 4, lty = 2)
legend(x = 2.3, y = 0.8, lwd = c(3, 3), col = c('gray', 'black'),
       legend = c('Admission Probability', 'Costs / Stakes'), 
       bty = 'n', cex = 1.3)
dev.off()